Coulomb stress accumulation along the San Andreas Fault system

Stress accumulation rates along the primary segments of the San Andreas Fault system are computed using a three-dimensional (3-D) elastic half-space model with realistic fault geometry. The model is developed in the Fourier domain by solving for the response of an elastic halfspace due to a point vector body force and analytically integrating the force from a locking depth to infinite depth. This approach is then applied to the San Andreas Fault system using published slip rates along 18 major fault strands of the fault zone. GPS-derived horizontal velocity measurements spanning then entire 1700 x 200 km region are then used to solve for apparent locking depth along each primary fault segment. This simple model fits remarkably well (2.43 mm/yr rms misfit), although some discrepancies occur in the Eastern California Shear Zone. The model also predicts vertical uplift and subsidence rates that are in agreement with independent geologic and geodetic estimates. In addition, shear and normal stress along the major fault strands are used to compute Coulomb stress accumulation rate. As a result, we find earthquake recurrence intervals along the San Andreas Fault system to be inversely proportional to Coulomb stress accumulation rate, in agreement with typical co-seismic stress drops of 1-10 MPa. This 3-D deformation model can ultimately be extended to include both time-dependent forcing and viscoelastic response. Smith and Sandwell Page 2 1/21/03 Introduction The San Andreas Fault (SAF) system, spanning over 1700 km from the Mendocino Triple Junction in the north, to the Gulf of California in the south, defines the complex tectonic boundary between the Pacific and North American plates. As the two plates slide past each other, the SAF system accommodates approximately 35-50 mm/yr of relative plate motion that is distributed across a 200 km wide zone [Working Group on California Earthquake Probabilities (WGCEP), 1995, 1999]. The SAF system is comprised of an intricate network of subfaults, each of varying geometry, locking depth, and associated failure properties. Earthquake recurrence intervals also vary dramatically along the SAF system subfaults, ranging from 20 years to over 300 years. In order to better understand the earthquake cycle and also help constrain faulting models of the San Andreas Fault system, geodetic measurements of interseismic, postseismic, and coseismic deformation are continually collected of the entire North American-Pacific plate boundary. While many previous studies of the SAF region have developed local fault-slip models to match regional geodetic observations of surface displacement [Savage and Burford, 1973; Savage et al., 1979; King et al., 1987; Li and Lim, 1988; Eberhart-Phillips et al., 1990; Savage, 1990; Lisowski et al., 1991; Feigl et al., 1993; Savage and Lisowski, 1993; Freymueller et al., 1999; Burgmann et al., 2000; Murray and Segall, 2001], our objectives are somewhat different in that we investigate the steady-state behavior of the entire San Andreas Fault system. By constraining relative plate motion, maintaining appropriate fault geometry, and implementing geodetic measurements spanning the entire system, we are able to model 3-D deformation and calculate stress accumulation. First, we investigate whether a single far-field plate velocity can be partitioned among parallel strands in order to accurately model near-field geodetic measurements. Second, we establish spatial variations in apparent locking depth along the main segments of the SAF system. Finally, we use our model to estimate secular buildup in Coulomb stress within the seismogenic layer and accumulation of scalar seismic moment. The primary purpose for developing our model is to estimate Coulomb stress accumulation rate and to explore its relevancy to earthquake occurrence and failure potential. Following the assumption that major earthquakes typically produce stress drops on the order of 1-10 MPa, estimates of Coulomb stress accumulation rate can provide an upper bound on the recurrence interval of a particular fault segment. Furthermore, recent studies of induced Coulomb stress changes propose that earthquakes may be triggered by stress changes as small as 0.1 MPa [King et al., 1994; Stein et al., 1994; Fialko and Simons, 2000; King and Cocco, 2001; Zeng, 2001]. High Coulomb stress accumulation rate has also been linked to areas of surface creep [Savage and Lisowski, 1993]. A better understanding of such stress-release processes at major plate boundaries, along with estimates of seismic moment magnitude, have also been the focus of recent earthquake hazard potential studies [WGCEP, 1995, 1999; WGNCEP, 1996]. In this analysis, we inspect the role of locking depth, fault geometry, and paralleling fault strands on accumulating interseismic stress along San Andreas Fault system, and investigate how such accumulation is related to shallow fault creep, earthquake recurrence interval, and seismic moment accumulation. Fourier Solution to 3-D Body Force Model For the last several decades, the most commonly used analytic models of fault-induced deformation have been based on the dislocation solutions of Chinnery [1961, 1963], Rybicki [1973], and Okada [1985, 1992]. The latter provide analytic expressions for stress, strain, and Smith and Sandwell Page 3 1/21/03 displacement in an elastic half-space due to a displacement discontinuity. While these dislocation models are accurate and computationally efficient when applied to individual faults or small fault systems, they may become computationally prohibitive when representing fault geometry over the entire North American-Pacific plate boundary. For example, 4x10 model calculations are required for 1000 GPS measurements and 400 fault patches. Modeling of InSAR observations could easily require 4x10 model calculations. However, if model calculations are performed in the spectral domain, the computational effort is substantially reduced. Rather than calculate the Fourier transform of the analytic solutions mentioned above, we instead solve the 3-D elasticity equations in the wave-number domain and then inverse Fourier transform to obtain space domain solutions. The key elements of our model derivation are summarized in Appendix A. In the two-dimensional case, our model matches the classical arctangent solution of Weertman [1964], both analytically and numerically (Appendix A). While this elastic half-space model currently ignores crustal heterogeneities and does not explicitly incorporate non-elastic rheology below the brittle-ductile transition, it produces reasonable estimates of first-order tectonic features comparable to other simple models [Savage and Burford, 1973]. To summarize our analytic approach (Appendix A), the elasticity equations are used to derive a set of transfer functions (in the wave-number domain) for the 3-D displacement of an elastic half-space due to an arbitrary distribution of vector body forces. The numerical components of this approach involve generating a grid of force couples that simulate complex fault geometry, taking the 2-D Fourier transform of the grid, multiplying by the appropriate transfer function, and finally inverse Fourier transforming. The force model must be designed to match the velocity difference across the plate boundary and have zero net force and zero net moment. There is a similar requirement in gravity modeling, where mass balance is achieved by imposing isostatic compensation and making the grid dimensions several times larger than the longest length scale in the system (e.g., flexural wavelength or lithospheric thickness). For this fault model, the characteristic length scale is the locking depth of the fault. The key is to construct a model representing a complicated fault system where the forces and moments are balanced. Our numerical approach is as follows: i) Develop a force couple segment from the analytic derivative of a Gaussian approximation to a line segment as described in Appendix A; this ensures exact force balance. For an accurate simulation, the half-width of the Gaussian must be greater than the grid size but less than the locking depth. ii) Construct a complicated force couple model using digitized fault segments. For each segment, the strength of the couple is proportional to the long-term slip rate on the fault segment and the direction of the couple is parallel to the overall plate boundary direction (not the local fault direction). This simulates the far-field plate tectonic force couple. Because the model has force couples, the vector sum of all of the forces in the model is zero but there is a large unbalanced moment because all of the force couples act in the same direction. iii) Double the grid size and place a mirror image of the force couple distribution in the mirror grid so the moment due to the image fault exactly balances the moment due to the real fault. Following these steps, we combine both analytic and numeric approaches to elastic fault modeling for analysis of the San Andreas Fault system. Modeling the San Andreas Fault Zone We apply our semi-analytic model (Appendix A) to the geometrically complex fault setting of the SAF system. After digitizing the major fault strands along the SAF system from geologic maps [Jennings, 1994] into over 400 elements, we group the elements into 18 fault segments Smith and Sandwell Page 4 1/21/03 spatially consistent with previous geologic and geodetic studies (Figure 1). Fault segments include the following regions: Imperial (1), Brawley (defined primarily by seismicity [Hill et al., 1975] rather than by mapped surface trace) (2), Coachella Valley-San Bernadino Mountains (3), Borrego (includes Superstition Hills and Coyote Creek regions) (4), Anza-San Jacinto (includes San Bernadino Valley) (5), Mojave (6), Carrizo (7), Cholame (8), Parkfield Transition (9), San Andreas creeping (10), Santa Cruz Mountains-San Andreas Peninsula (11), San Andreas North Coast (12), South-Central Calaveras (13), North Calaveras-Concord (14), Green Valley-Bartlett Springs (15), Hayward (16), Rodgers Creek (17), Maacama (18). The fault system is rotated about its pole of deformation † (52 N, 287 ° W) into a new co-ordinate system [Wdowinski et al., 2001], and fault segments are embedded in a 1-km grid of 2048 elements along the SAF system and 1024 elements across the system (2048 across including the image). We assume that the system is loaded by stresses extending far from the locked portion of the fault and that locking depth and slip rate remain constant along each fault segment. Each of the 18 SAF segments is assigned a deep slip rate based on geodetic measurements, geologic offsets, and plate reconstructions [WGCEP, 1995, 1999]. In some cases, slip rates (Table 1) were adjusted (+/5mm/yr on average) in order to satisfy an assumed far-field plate velocity of 40 mm/yr. This constant rate simplifies the model and, as we show below, it has little impact on the near-field velocity, strain-rate, and Coulomb stress accumulation rate. Moreover, it provides a remarkably good fit to the geodetic data, except in the Eastern California Shear Zone, where misfit is expected due to omission of faults in this area. Because slip estimates remain uncertain for the Maacama and Bartlett Springs segments, we assume that these segments slip at the same rates as their southern extensions, the Rodgers Creek and Green Valley faults, respectively. After assigning these a priori deep slip rates, we estimate lower locking depth for each of the 18 fault segments using a least squares fit to 1099 GPS-derived horizontal velocities. Geodetic data for the southern SAF region, acquired between 1970-1997, were provided by the Crustal Deformation Working Group of the Southern California Earthquake Center (SCEC Velocity Field, Version 3.0). GPS velocities for the Calaveras-Hayward region were provided by the U.S. Geological Survey, Stanford University, and the University of California, Berkeley and reflect two data sets, one of campaign measurements (1993-1999) and one of BARD network continuous measurements. Data used to model the northern region of the SAF system were obtained from Freymeuller et al. [1999], and represent campaign measurements from 1991-1995. These four geodetic data sets combine to a total of 1099 horizontal velocity vectors spanning the entire San Andreas Fault zone (Figure 2a). Geodetic Inversion The relationship between surface velocity and locking depth is nonlinear, thus we estimate the unknown depths of the 18 locked fault segments using an iterative, least-squares approach based on the Gauss-Newton method. We solve the system of equations Vgps(x,y) =Vm(x,y,d), where † Vgps is the geodetic velocity measurement, † Vm is the modeled velocity, and d is the set of locking depth parameters that minimize the weighted residual misfit † c 2 . The data misfit is Smith and Sandwell Page 5 1/21/03 † Vres i = Vgps i -Vm i